Predicting the splash of a droplet impinging on solid substrates

The impingement behaviours of droplets towards solid substrates depend on the liquid properties, impingement velocity and solid surface conditions, such as wettability and roughness. However, the prediction regarding whether the droplet splashes after the impingement, is still an open question. Here we show that the splashing can be predicted by the pressure balance of the liquid film appearing beneath the impingement droplet coupled with the modified energy balance equation. Hydrodynamic and hydrostatic pressures are the driving forces for the droplet’s radial spreading, while the capillary pressure at the rim edge and viscous stress oppose the driving forces. Thus, splashing occurs when the driving forces overcome the opposing forces. Moreover, the splashing condition is affected by various surface factors, such as wettability and surface roughness. Our work would pave the way to understand the basic physics for rim or liquid film fragmentation and enabling advances in important for engineering field such as printing, sprays for cooling and pesticide.

www.nature.com/scientificreports/ liquid film behaviour has helped increase our understanding of the mechanism behind the splashing behaviour of the droplets, the researchers have still not reached a consensus in their findings [21][22][23][24][25] .
Recently, Yonemoto and Kunugi 26 proposed an analytical model for the prediction of the maximum spreading contact-area diameter of the droplets based on the energy balance approach. The model was in good agreement with the experimental data covering various liquid viscosities and surface tensions 27 . This model suggested that the droplet ideally spreads over smooth solid surfaces without splashing. However, for practicality, the understanding of the droplet wetting behaviour on rough solid surfaces, including deposition and splashing, is very important [18][19][20][26][27][28] . The spreading and retracting processes of the droplet mainly depend on its inertia and wettability, complicating the evaluation of the maximum spreading area 29 . For both smooth and rough surfaces, droplet splashing occurs during its spreading. Therefore, in the present study, we have modified the energy balance equation for both smooth and rough solid surfaces and developed a novel model for predicting the splashing condition of a droplet. This model combines the pressure balance of the liquid film with a simplified circumferential-instability model. Finally, the splashing criterion is obtained by solving the modified energy balance equation combined with the newly developed pressure balance relation. The developed model can predict the splashing condition on both smooth and rough solid substrates. In addition, the modified energy balance equation can also predict the spreading contact-area diameter in the deposition region on both the substrates.

Methods
Water-ethanol binary mixtures were used in our study. The mass concentration was varied from 0 to 99.4 wt.% ethanol (pure ethanol; Nacalai Tesque, Inc., 99.4 wt.%). Here, the viscosity of the water-ethanol binary mixtures exhibits the maximum value at around 40 wt.%. Therefore, we chose the weight percentage of 6 cases to investigate the effect of the surface tension on the splashing behaviour under similar liquid viscosity µ l . Concretely speaking, the ultrapure water (Wako Pure Chemical Industries, Ltd., σ lg = 0.0719 N m − 30 . Polycarbonate (PC) was used as the solid material. To evaluate the effect of the surface roughness on the droplet splashing behaviour, the solid material was polished with a grinder-polisher (MetaServ TM 250 Grinder-Polisher, Buehler Ltd., Lake Bluff, IL, USA). The surface was prepared by using three types of polishing sheets, 408-400AU (Grit #400), 408-240AU (Grit #240) and 408-120AU (Grit #120), Sankei Co. Ltd., JAPAN. Silicone rubber (SR) substrate was also used as an additional hydrophobic substrate with negligible adsorption. The surface morphological conditions were measured using a laser scanning microscope (LEXT OLS4100, Olympus Co. Ltd., JAPAN). The arithmetical mean-roughness value (R a ) was measured for each substrate as 0.03 µm for the bare plate, 0.33 µm for #400 plate, 0.99 µm for #240 plate and 1.25 µm for #120 plate. R a for SR was 0.02 µm. Prior to the experiment, the solid materials were rinsed with ethanol and purified water for the SR substrate and the purified water for the PC substrates. Then, the materials were dried with an air blower.
The droplets were gently released without any initial velocity using a micro-syringe from different heights (5-2400 mm). The droplet volume ranged from 3.4 to 6.6 µL. The errors in the droplet volumes of the droplets were within 3% in the present study. In addition, the difference between the measured vertical and horizontal initial droplet diameters were within 10% 31 , indicating a good repeatability of the experiment. A high-speed video camera (HX-5, NAC image technology, Ltd., Japan) with a microscope (Leica Microsystems, Welzlar, Germany) or a micro-lens (Nikon AF-S VR Micro-Nikkor 105 mm f/2.8G IF-ED) was used to capture the impingement behaviour with the frame rate of 20,000 fps. Each condition was repeated three times. The critical splashing velocity of a droplet was determined by observing the secondary droplet falling from the same height thrice. The validity of the present experiment was confirmed by comparing with the existing splashing models (Section S5 in the supplementary information). Notably, the identification of the splashing type, prompt or corona, is out of scope for the present study [32][33][34][35] . The temperature and relative humidity were maintained within 20.0-25.0 ℃ and 50.0-55.0%, respectively.

Experimental results
Spreading of droplets on smooth and rough solid surfaces. In the droplet impingement process, the droplet shape exhibits a damped vibration process as the time advances. The droplet shape change also affects the behaviour of the contact-area diameter d cont (t), here, t is the time. For weak hydrophilicity, the change in the contact-area diameter post-impingement reaches 0, which corresponds to the stationary point represented by dd cont (t stp )/dt ≈ 0, here, t stp is the time at the stationary point. Hereafter, the contact-area diameter d cont (t stp ) or radius r cont (t stp ) at the stationary point is represented by d stp or r stp . The size of the contact-area diameter at the second stationary point is usually smaller than the first one. However, for strong hydrophilicity, because of the intermolecular interaction or the effect of surface roughness on wettability, the behaviour of the contact-area diameter becomes complicated. Figure 1 shows the time evolution of the normalised contact-area diameter β(t) (= d cont (t)/d 0 : d 0 is the initial droplet diameter) for (a) 0 wt.% and (b) 99.4 wt.% on the bare plate and (c) 0 wt.% and (d) 99.4 wt.% on #400 plate. The droplet is released from different heights of 5, 10, 50 and 100 mm. Figure 2 depicts the earlier stage in Fig. 1 and shows the time evolution of β(t) post droplet impingement. In Fig. 2, the red points represent when the change in β(t) approximately becomes 0 (i.e. the first stationary point of dβ(t stp )/ dt ≈ 0). Here, the β(t stp ) is the normalised contact-area diameter at the stationary point. Hereafter, β(t stp ) is represented by β stp . Evidently, for (a), the contact-area diameter reaches the stationary point just after the droplet impingement for each release height. Subsequently, β(t) becomes small (as in the retraction process) and reaches www.nature.com/scientificreports/ a constant value. For each release height, the stationary point assumes the maximum value because β(t) after the first stationary point does not exceed the value of β(t stp ). In the case of (b), β(t) reaches the stationary point just after the impingement and exhibits the retraction process for each release height (Fig. 2). However, after retracting, β(t) gradually increases again and becomes larger than the value at the first stationary point. Basically, the static wetting condition of 99.4 wt.% ethanol exhibits complete wetting of the bare plate 36 . This means that the droplet continues to spread over the solid surface during the droplet impingement process. Therefore, the adsorption effect on the contact line motion becomes prominent after the kinetic energy is almost consumed. This result indicates that there are two processes for the droplet spreading 29 . One involves inertia wherein the spreading is promoted by the fluid motion. The other one is wettability dominated process in which the spreading is promoted by the molecular scale interaction. Moreover, the surface morphology also affects the spreading, as shown in Figs. 1c, d and 2c, d. In these cases, the droplet behaviour becomes more complicated than that in the bare substrate because the surface roughness induces the pinning of the contact line, making it hard to recede. In the (c) and (d) cases, the retraction of the contact line after the first stationary point becomes small. Therefore, β(t) after the first stationary point exhibits various behaviours depending on the release height, wettability and surface roughness. In other words, it depends on the behaviour of the contact-area diameter after the first stationary point whether the stationary point becomes the global maximum, local maximum, or inflection point. If the stationary point contains the largest value through the whole droplet impingement process, the spreading contact-area diameter becomes the maximum spreading contact-area diameter. At least, the droplet spreading behaviour until reaching the first stationary point is the inertia dominated process. Therefore, there is no physical contradiction to the prediction of the spreading contact-area diameter of the droplet at the first stationary point using the energy balance approach (Section S4). Figure 3 shows the splashing images of water-ethanol binary-mixture droplets on the bare (smooth solid surface) and #240 PC substrates. For bare substrate, the splashing velocity of the water droplet is larger than that of other liquids. Conversely, for #240 PC substrate, the differences in the splashing velocities are small in terms of the liquid concentration. For low concentrations (0 wt.% and 5 wt.%), the splashing velocities on the bare substrate are larger than those on other roughened solid substrates. Similar tendency has been observed for water in another study 20   www.nature.com/scientificreports/

The model for splashing
Pressure balance of the liquid film. The typical image for the droplet impingement process reflects that the liquid film spreads over the solid surface just after the droplet impingement, and the tip of the liquid film reaches the first stationary condition dd cont (t stp )/dt ≈ 0, where the contact-area diameter d cont (t) achieves the maximum value and the droplet exhibits retraction. Finally, the droplet reaches equilibrium through the damped vibration process. The droplet splashing occurs in a series of the impingement process if the kinetic energy is sufficient to eject the secondary droplets from the tip of the liquid film. As shown in Fig. 4a, the existence of the liquid film when the droplet impinges on the solid substrate is very important for the discussion of the splashing behaviour 17 . Especially, it can be assumed that the splashing occurs when the secondary droplets are generated from the tip of the disc-like liquid film. The pressure balance is considered for this concept 37 . At a moment in the droplet impingement, the liquid starts spreading in the radial direction, and the hydrodynamic pressure, P hyd , driving the liquid film in the radial direction is expressed using the initial radial velocity u H 0 as follows: In Eq. (1), ρ l is the liquid density and the initial radial velocity u H 0 is assumed to be 3u/8, where u is the impinging velocity perpendicular to the solid substrate. u H 0 (= 3u/8) can be estimated by considering mass conservation where a cylindrically shaped droplet of diameter d 0 is assumed for the sake of convenience 26 . Then, the hydrostatic pressure P stat , another driving force is expressed as follows: In this equation, d 0 represents the initial droplet diameter and g is the gravity. Then, if a thin disc-like liquid film protrudes in the radial direction at a moment in the impingement, the surface tension exerted on the tip of the liquid film (f tip ), as shown in Fig. 4a, b, is expressed as f tip = (2πr film σ lg /2) + (2πr 0 σ lg /2) characterized by r film and r 0 (= d 0 /2), where σ lg is the surface tension between the liquid and gas. Then, a simplified model of the circumferential instability to the liquid film is applied in the study, as shown in Fig. 4c, where the concave (L conc ) and convex (L conv ) arc lengths are approximated the same: L conc ≈ L conv (Section S1 in the supplementary information). A secondary droplet is assumed to eject from every convex position along the circumference based on the simplified circumferential-instability model. From this assumption, the effective area A eff on the tip of the (2) P stat = ρ l gd 0 . www.nature.com/scientificreports/ liquid film where the surface tension exerts can be modelled as h film × (2πr 0 × 1/2). Finally, the Laplace pressure P Lap exerted on the tip of the liquid film can be derived from f tip /A eff as follows: In addition, the viscous stress P vis is expressed by    (4), the secondary droplet ejected from the tip of the liquid film is generated when the following relation holds: The boundary condition whether the droplet splashed or not is P hyd + P stat = P Lap + P vis . A precise prediction of the size of the secondary droplet ejected from the liquid film is very difficult. Therefore, in the present study, the secondary droplet size is characterised by the droplet size of h stp evaluated by postulating a disc-shaped droplet as V 0 /πr stp 2 , where r stp is the spreading contact-area radius at the first stationary condition, V 0 is the volume of droplet and h stp is the droplet height at r stp 26 . Here, the first stationary condition represents the situation in which the droplet spreads just after the impingement and the change in the contact-area diameter d cont (t) approximately becomes zero (dd cont (t stp )/dt ≈ 0). The droplet spreading behaviour until reaching the first stationary point is the inertia dominated process ("Spreading of droplets on smooth and rough solid surfaces" section). Thus, if the contact-area diameter in the first stationary condition holds the largest value throughout the droplet impingement process, the spreading contact-area diameter is called the maximum spreading contact-area diameter d m . It depends on the behaviour of d cont (t) after the first stationary condition whether the stationary point becomes the global maximum, local maximum or inflection point. Furthermore, the secondary droplet size is affected by the wettability because the spontaneous capillary-driven spreading will occur when the droplet contacts the solid substrate 38,39 . Even in the case of an impact of a solid body in a liquid bath, the wettability affects the splashing behaviour of the liquid 40 . Therefore, based on the above concept, the thickness of the liquid film h film can be modelled as follows (Section S2 in supplementary information): (5) P hyd + P stat ≥ P Lap + P vis . Effect of surface roughness on spreading. Wettability of droplet on the solid substrate is affected by the surface morphology and the interaction between the solid and liquid molecules. The intermolecular interaction is an essence of the physics of wettability. Therefore, the adsorption of liquid molecules on a solid substrate is an important phenomenon 43 . The difference in the strength of intermolecular interaction would be indirectly exhibited by an indicator for wettability, such as contact angle. For example, the adsorption of liquid molecules on a hydrophobic substrate 44 is not as significant as that on a hydrophilic one 36 . This is inevitable for the energy balance concept used for the prediction of the spreading factor. Thus, in the present study, the conventional energy balance equation is also modified based on the same concept.

Ejection of secondary droplets
Morphological effect on radial spreading velocity. As the surface roughness increases, the effect of the solid surface on the fluid motion would also increase, inducing changes in the viscous force and dissipation. This type of morphological effect must be considered in the pressure balance model and conventional energy balance equation 26 . The evaluation of the roughness effect strictly on the fluid motion was difficult although the relationship between the radial spreading velocity and surface conditions, such as wettability and roughness, was highlighted 20 . Therefore, in the present study, the morphological effect is evaluated by considering the relationship between the work of adhesion for solid-liquid and the hydrodynamic force radially pushing the droplet. Therefore, hydrodynamic force comparable to the work of adhesion is at least needed to push the liquid forward in a radial direction. For example, for a bare substrate, the relation among the initial spread diameter d bare , the initial radial velocity u H 0,bare and the contact angle θ bare can be evaluated as follows: For a rough substrate, the relation among the initial spread diameter of d rough , initial radial velocity of u H 0,rough and contact angle of θ rough can be considered as follows: By taking the ratio of Eqs. (7)-(8), the following relation is deduced: Equation (9) implicitly indicates the morphological effect on the fluid motion. Therefore, if u H 0,bare /u H 0,rough is defined as k, it is related to the radial spreading velocity in Eqs. (1) and (4). Thus, from Eqs. (1)-(5), the following equation can be derived: where β stp is the factor of the spreading contact-area diameter at the first stationary condition. In the present model, We and Re are evaluated by the initial droplet diameter d 0 like as We = ρ l u 2 d 0 /σ lg and Re = ρ l ud 0 /µ l , respectively. Similarly, k is incorporated in the viscous dissipation term in the energy balance equation as follows: where µ l is the liquid viscosity and t stp is the time at which r stp is reached. Thus, the viscous dissipation term in the conventional energy balance equation 26 is replaced by Eq. (11).
Adsorption effect on spreading. In the present study, polycarbonate and water-ethanol binary-mixture liquids were used as the solid material and the test fluid, respectively. In this combination, for low-ethanol concentration, the more the surface roughness of the solid substrate, the less is the wettability. In contrast, if the ethanol concentration is high, the more the surface roughness of the solid substrate, the more is the wettability of the liquid 36 . However, even for bare substrate, more the ethanol concentration, the more is the wettability. In addition, the contact line would be easy to move in the advancing direction if the wettability exhibits hydrophilicity 45 . These phenomena indicate that the deformation of liquid surface can contribute to the contact line movement (radial spreading). Then, the additional work infiltrating into the surface morphology, such as protrusions or (6)   1 + cos θ bare 1 + cos θ rough . www.nature.com/scientificreports/ grooves would generate. In particular, polycarbonate has a hydrophilic tendency for water-ethanol binary-mixture liquids, making the adsorption important for the wetting behaviour. In microscopic point of view, the liquid molecules in the vicinity of the solid surface would behaves like a solid [46][47][48] because the bonding energy of solidliquid molecules will be larger than that of liquid-liquid molecules in the case of the PC substrates, although it would depend on the combination of solid and liquid in an actual situation and would be also influenced by the surface morphology. The effect of the solid-liquid adsorption on the liquid molecules motion will gradually decreases with increasing distance from the solid surface 49 . Thus, it can be assumed that the adsorption at the solid-liquid interface would be categorized into two parts. One is the strong adsorption effect on the liquid molecules in the vicinity of the solid surface. The other is the weak adsorption effect on the liquid molecules in the bulk side near the solid surface. The latter can consider as the adhesion energy characterised by the macroscopic information such as the contact angle and the surface tension. The former will be more microscopic interaction than the latter and cannot express by the macroscopic information. However, it would be particularly important in the dynamic process of the wetting behaviour such as the slip condition 50 and the dynamic contact angle [51][52][53][54][55] . Therefore, the energy terms related to the adsorption (E ads ) and the infiltration into the surface morphology (E infil ) are added into the energy balance equation 26 .
where E kine , E surf , E grav , E sprd and E def are the kinetic energy, initial surface energy, gravitational potential of the droplet, adhesion energy and deformation energy after the impingement, respectively, and the additional energy E additional is considered as follows: Here, it is assumed that the energy can be indirectly estimated from E def because the adsorption and the surface morphology affect the droplet surface deformation through the contact line motion. The concrete expression for E ads is modelled as follows: In Eq. (14), the ratio of E ads to E def , affected by the intermolecular interaction between the solid and liquid, is represented as g ads . As to g ads , Langmuir-type function is postulated using the surface tension of liquid (σ lg ) and the critical surface tension of solid (σ c ) based on the liquid molecule adsorption to the solid substrate as follows 36 : where A and B are arbitrary parameters. In Eq. (14), the negative value for g ads means that the reduction of energy for the surface deformation corresponds to the relative increase of other energies in the energy balance equation, which leads to the spreading of droplet. The positive value indicates that the energy is consumed in the surface deformation rather than the spreading. Here, the critical surface tensions σ c for each PC substrate are 0.0247 N m −1 for bare, 0.0276 N m −1 for #400, 0.0303 N m −1 for #240 and 0.0317 N m −1 for #120, respectively 36 . The value of σ c for SR substrate is 0.0196 N m −144 .
As for E infil , the following relation is assumed, In Eq. (16), g infil is defined by the following relation as the synergy effect of the intermolecular interaction and the surface morphology on the wettability: In this equation, C and D are arbitrary parameters. f represents the relative surface-roughness area 36 . Here, the values f for PC substrates are 1 for bare plate, 1.12 for #400, 1.23 for #240 and 1.28 for #120, respectively. In the case of the SR substrate, the value of f is treated as 1. E infil reflects the work for the infiltration of liquid into the surface morphology, and g infil (= E infil /E def ) qualitatively expresses the wetting states, such as Wenzel or Cassie state (Section S3 in supplementary information). The parameters A-D in Eqs. (15) and (17) Figure 5a shows that the critical We becomes small as the relative roughness R a /d 0 increases. This tendency is related to the static wetting condition 36 . As the surface roughness increases, the hydrophobicity of water droplet (0 wt.%) also increases, corresponding to the Cassie situation in which an air pocket exists between the liquid and solid. This tendency also appears in the spreading behaviour at the stationary condition (Fig. S3). Therefore, the incomplete wetting situation would induce a morphological effect on the circumferential instability of the liquid film and results in the low critical We for splashing. Considered the static wetting condition, the case for (b) 5 wt.% can be also considered to be following the same mechanism. Conversely, for moderate ethanol concentrations of (d) 40 wt.% and (e) 70 wt.%, the critical We increases as the relative roughness increases. In these cases, the wettability of the static droplet becomes large as the surface roughness increases 36 . Therefore, the increase in the adhesion to the solid substrates results in the increase of the critical We for the splashing. The case for (c) 20 wt.% would be a transition region from (a) to (d  26,27 . These results prove the validity of the present model. Figure 7 shows the liquid concentration dependence on the ratio of non-dimensional energies for the capillary (R cap ) and the viscous dissipation (R vis ) for each solid substrate under the splashing condition. R cap and R vis are E cap /E ini and E vis /E ini , respectively. Here, E cap = E sprd + E def + E additonal and E ini = E kine + E surf + E grav 27 . For the bare substrate (Fig. 7a), the splashing occurs in the viscous dissipation dominated region for all liquid concentrations. However, in the low-ethanol concentration region, R cap gradually increases with the surface roughness and splashing occurs in the capillary dominated region for the #120 substrate for water and 5 wt.% ethanol. These results  www.nature.com/scientificreports/ indicate that the critical We is small for the splashing in the capillary region. Conversely, the effect of wettability on the splashing behaviour becomes small as the hydrophilicity increases (high-ethanol concentration region). The difference in the critical We becomes small as the ethanol concentration increases (Fig. 5). Similar behaviour would be observed for metal substrates because the surface free energy of the high-surface energy solid is larger than that of the low-surface energy solid, such as polymer or plastic substrate. In the present study, the maximum and minimum values of the critical capillary numbers (Ca) where the splashing occurs are 0.23 and 0.03, respectively. Therefore, the effect of wettability on the splashing behaviour may disappear if the critical Ca exceeds 0.23 in the present study 56 . Notably, this tendency may depend on the combination of the solid surface morphology/wettability and liquid property/composition.

Conclusion
In the present study, the droplet splashing behaviours for water-ethanol binary-mixture liquids to some lowsurface-energy solids were experimentally investigated. Then, the theoretical model for the prediction of the splashing conditions was developed, considering the pressure balance for the liquid film. In addition, the energy balance equation was modified for both smooth and rough solid surfaces. The modified energy balance equation can predict the spreading factor on the rough solid surface in the deposition region. Furthermore, analytical results obtained by solving the pressure balance equation combined with the modified energy balance equation are in good agreement with the experimental data for the critical Weber number for splashing. The result indicates that the splashing conditions depend on the wettability between the solid and liquid, and the solid surface roughness in addition to the viscosity. In particular, the effect of adsorption on the wettability is an important factor. The present model mainly focuses on the processed surface that is not so rugged, and relatively smooth surface (processed by sandpaper), which corresponds to We ε ≤ 1 based on the classification by Garcia-Geijo et al. 20 . However, prior to our study, no theoretical model could predict both the splashing conditions for the droplets and the spreading contact-area diameter at the first stationary condition, such as the maximum spreading diameter on the smooth and rough solid surfaces. Although the present model cannot apply to the rugged solid surface where the contact angle is hard to define, predicting the splashing condition on such rugged solid surfaces could be useful as the basic model or methodology.